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Abstract 

We have done a finite-size scaling study of a continuous phase transition altered by the quenched 



bond disorder, investigating systems at quasicritical temperatures of each disorder realization by 

using the equilibriumlike invaded cluster algorithm. Our results indicate that in order to access the 

c- 1 

thermal critical exponent y T , it is necessary to average the free energy at quasicritical temperatures 

of each disorder configuration. Despite the thermal fluctuations on the scale of the system at 

the transition point, we find that spatial inhomogeneities form in the system and become more 

pronounced as the size of the system increases. This leads to different exponents describing rescaling 

f—i ' of the fluctuations of observables in disorder and thermodynamic ensembles. 
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The question of how disorder modifies a phase transition is still not fully resolved, notwith- 
standing decades of research. Many issues remain open, even in the case of the ferromagnetic 
transition in the presence of quenched disorder, from the qualitative nature of the ordered 
phase near criticality |l|, |2j, the presence of Griffiths singularities |3j, the role of multiple 
length scales [4J, to violation 5| of the proposed 6] bound on the thermal critical exponent. 

An important aspect present in the majority of the cited problems is the non-self- 
averaging effect between different disorder realizations [7j, leading to the question of proper 
averaging over disorder. The standard approach is to reduce the problem to an effective 
translationally invariant one, by averaging the free energy over disorder configurations at 
a given temperature [8|. Such averaging permits one to apply the perturbative renormal- 
ization group (RG), and crucially simplifies the finite-size scaling (FSS) studies based on 
numerical simulations as well. When applied to a classical ferromagnetic model with bond 
randomness, displaying a continuous phase transition in the pure case such as the one we 
consider in this paper, the conclusion seems to be [9| that the disorder, if relevant in the 
sense of the Harris criterion 10] , merely changes the values of the critical exponents, while 
the nature of the transition remains similar to that of the pure case. 

We challenge this notion with a FSS study in which the average is taken at quasicritical 
(i.e. finite size critical) temperatures of each disorder sample. Except for some partial 
attempts, a systematic study using this type of averaging is still missing, since it requires an 



except iona 
algorithm 



numerical effort. The recently proposed equilibriumlike invaded cluster (EIC) 



11 



12J enables such a procedure to be technically feasible. 



There exist prior works revealing conceptual problems in the transition with bond dis- 



of observables at criticality 
141 . Correlation length has 



order. It has been pointed out that the self- averaging ratios [7 

depend on averaging 13| although they should be universal 

also been shown to be non-self-averaging near criticality 15|. Pazmandi et. al. |5| have 

shown that standard averaging obscures the intrinsic thermal exponent y T = 1/V, when it is 

superior to the exponent y T describing disorder fluctuations of quasicritical temperatures, 



Vr < Vt- 



(1) 



In such a case, in standard averaging, y~ T always takes over v T i n FSS relations, and the 
information about y T is lost. It has therefore been suggested [13, lla] that in order to find 
the true critical exponents one has to perform a disorder average of observables, taken at 



the quasicritical temperature of each disorder sample. 



We use the E 
bond disorder 



C approach on the two-dimensional (2D), three-state Potts model with 



171 ]. where disorder is relevant by the Harris criterion, since the specific-heat 



critical exponent of the pure model is positive (a = 1/3). Our main result shows inequality 
([T]) by an explicit numerical calculation, justifying the assumption by Pazmandi et al. [5] in 
a classical system. We suggest an alternative interpretation of the lack of self-averaging as 
an emergence of frozen inhomogeneities in observables at all scales, which may be analyzed 
further from the statistics of the largest cluster at criticality. 

Since the Chayes bound [6j, i.e. y T < d/2, is exclusively a result of averaging at a 
unique temperature, it is irrelevant for the intrinsic thermal exponent y T . For this reason 
our findings might a 



was observed 



18 



so be of interest to systems in which a violation of the Chayes criterion 
jN20l|. Experimentally, the intrinsic exponent y T would be observed only in 
direct calorimetric measurements of the heat capacity, while the indirect measurements (e.g. 
birefringence 21] ) would yield the exponent tjj- obeying the Chayes criterion. 



We introduce disorder in the Potts model [22j as a random bond dilution 



H = E-^OW,- 1 ); 



(2) 



<i,j> 



where random couplings Jy are restricted to neighboring lattice sites i, j and take zero 
value with concentration c, or Jy = J otherwise, a denotes the Potts variable with q = 3 
discrete states. For the purpose of numerical simulations, the partition function of the Potts 
model is written in terms of the random cluster model, by using the Fortuin and Kasteleyn 
(FK) graph expansion [23] 



Z= ^p^).(l-p)^(7). g c(7) ; 



7er c 



where the bond probability 



V 



(3) 



(4) 



takes the role of temperature, while 6(7) and 0(7) denote the number of bonds and connected 
components "FK clusters" in the graph 7, respectively. The summation runs over the set of 
all possible FK graphs T a , compatible with the given disorder configuration a. 

When applied to a system with disorder, the EIC algorithm 11[ simulates it at the 
quasicritical (i.e., finite-size critical) bond probabilities p c a (L) (related to temperature by Eq. 



(jlj)) belonging to each disorder configuration a, defined by the onset of the percolation of 
the largest FK cluster. This property allows us to consider separately the scaling of thermal 
fluctuations of an observable around its thermodynamic mean and the disorder ensemble 
fluctuations of the thermodynamic means. Our algorithm is an extension of the invaded 



cluster (IC) algorithm 24J and has a similar mechanism of self-regulating to the quasicritical 



joint, but with the crucial difference of generating the equilibrium thermodynamic ensemble 



12] . The duration of a Monte Carlo (MC) step of the EIC algorithm is approximately 
the same as that of the IC algorithm, but because of the canonical constraint, the EIC 
algorithm requires a certain number of thermalization steps (never exceeding 5000 in the 
present paper). The EIC algorithm generalizes to the problem with bond dilution in a 
straightforward way since disorder merely excludes some configurations from the set of all 
possible FK graphs of the pure case. 

Throughout this paper we denote by [■] the disorder average taken at p c a (L) of each 
disorder configuration and by ~ thermodynamic average for a given a. Numerical results 
presented are based on simulations on square lattices of linear size L ranging from 64 up to 
896, with two disorder concentrations, c = 0.125 and 0.25, using the statistics of 400 and 600 
disorder configurations respectively, with 20000 Monte Carlo steps (MCS) per disorder con- 
figuration after thermalization. Calculation of a running average reveals that such a disorder 
statistics is sufficient to determine each [p^,(L)] to seven significant digits. The disorder is 
introduced microcanonically, i.e. the exact number of bond vacancies corresponding to a 
concentration is randomly distributed on the lattice. 

The thermal critical exponent y T is calculated from the magnetization-energy cumulant 



Tj-me _ m cx. e -a m a ' e a /r\ 

where m a and e a denote the order parameter and energy density of a disorder configuration 
ex, respectively. The values of y T in Table fl] are obtained by fitting (Fig. [T] a) the averaged 
data to the power law form [U™ e ] oc L Vr ~ d . The simple power law describes the scaling of 
[U™ e ] for c = 0.25 in the entire range of lattice sizes and for c = 0.125 only the data from 



the smallest size deviate [25|] . Since the scaling corrections are always found to be important 
in approaches using the averaging at a unique temperature (see e.g. [26]), their negligibility 
can only be attributed to the averaging procedure we used in this work. Exponent y T also 
shows a negligible concentration dependence, contrary to the previous studies where the 
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Figure 1: (color online) Slopes in (a)-(c) correspond to y T — |, y T — | and i — x respectively, pointing 
out the difference between the three exponents. The disorder statistics determines the precision of [C/™ e ] 
for each L to the fourth significant digit and the precision of 8p c and 8e c to the third significant digit. The 
error bars for L d / 2 [U™ e ] are smaller than the symbols. 



standard averaging was used [27 1. 



Table I: Rescaling exponents of [U™ e ], 8p c a and 8e c a . Each set of data has been fitted in the entire range of 
sizes to a single power-law except for the \U™ e ) for c = 0.125, where the smallest lattice has been excluded. 
Error bars have been estimated by jackknife binning. 




Vt 



d — x 



0.88(4) 1.35(5) 
0.95(2) 1.14(4) 



To understand how the exponent y T becomes obscured when the averaging at a unique 
temperature is used, we examine the fluctuations of quasicritical bond probabilities p c a (Fig. 
[T] b) in disorder ensemble (ensemble of random samples) 



Sf = VK 2 ] - [P c a ? oc L~ 



(6) 



By virtue of Eq. (J3J) p c obeys the same scaling law as the fluctuations of quasicritical 
temperatures T£. The fact that beyond statistical errors y T is dominant over y T (Tab. HI), 
leads to the conclusion that, by averaging at a unique temperature one effectively measures 
the exponent y T [5| . The range of values for different disorder concentrations that we have 
obtained for y T corresponds to the findings for the thermal exponent in previous studies, 



which used averaging at a unique temperature 28M30J for the same system. For example 



Jacobsen and Cardy 28| have estimated the value to 0.96(4). Unlike y T , the exponent y T 
displays a strong dependence on disorder, similar to that found in previous studies. 

From U™ e we have determined the singular part of the energy density which scales 
oc L Vr ~ d . Considering disorder ensemble fluctuations of thermodynamic means of energy 
density e£, taken at T£ (Fig. [T]c), 



5e c = ^/[e%}-K}^L-^ (7) 

we obtain that L d 5e c increases, with the exponent d — x (Tab. [T|) significantly larger than y T 
or y T for both c. Since Se c is the sum of the analytic and the singular part, the exponent x 
las to be attributed to the analytic part, which is dominant over the singular one. In Ref. 



281 ] the analysis of the energy fluctuations by the standard disorder averaging produces the 
exponent <i — 0.91(1) = 1.09(1), significantly larger than their value for the thermal exponent 
0.96(4). 

If the disorder had merely changed the values of the critical exponents found for the pure 
case, y T and d — x would have to be equal to y T , as determined by the singular part of the 
free energy. However, we show the exponents to be different (Fig. [p. 

The evidence presented leads to the conclusion that changing a disorder configuration a 
while keeping the concentration of disorder constant, drives the system out of the critical 
area. A consequence of this scenario is the formation of spatial inhomogeneities at all scales. 
Consider the free energy density for a single disorder configuration a 

/« = L- d f a {L yT (T - T^(L)), W»h) + T ■ e a (T, h; L), (8) 

where h is the magnetic field and Te a is the analytical part. Since we have demonstrated 
that inequality ([I]) applies, changing a produces a shift in the argument of expression (jHJ) 
oc L Vt ~^ t , which diverges for L — > oo. However, when a system of size L can be driven out 
of criticality by changing a, this shift is even larger for any of its parts of size V < L since 
jj-v-r > jj-vt Thg shift may be positive or negative in any of the parts since the disorder 
configurations in them are uncorrelated, so the system locally prefers either the low- or the 
high- temperature phase, resulting in spatial inhomogeneities. The contributions from local 
parts accumulate in the analytic part of the energy, causing the disorder ensemble fluctua- 
tions (i.e., the exponent x). If the opposite to inequality §\§ were true, the inhomogeneities 
would have disappeared in the limit L — > oo, since L Vt ~^ t would be at most of order 1. 
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Figure 2: (Color online) Spatial dependence of rn a (r) for a single disorder configuration at p% obtained 
with 10 6 MCS for L — 2048. The spatial average of rn a (r) in this disorder configuration is 0.29164 and is 
labeled by the tick on the scale. 

We demonstrate the existence of the spatial inhomogeneities on the example of the local 
order parameter m a (r), defined as the thermodynamic probability that the spin on the site 
r belongs to the largest FK cluster 

= E { .„ t> ^(q-^(W.W) 

Em,{6} w «({ct},{6}) 

where n a is 1 if <j(r) belongs to the largest FK cluster and otherwise. By w a we denote 
the statistical weight of a configuration in the joint space of spin and bond variables. The 



spatial average of m a (r) is equivalent to the standard Potts order parameter 31(]- In Fig. [2] 
rn a (r) in a single disorder sample with 25% of disorder at p c a (L) is presented. The behavior 
of fn a (r) contrasts with the situation in the pure model, where all the lattice sites are 
equivalent and m a (r) is constant. It illustrates how a fractal volume, pinned by disorder, 
emerges as the preferred ordering volume on which the local order parameter is close to 1. 

We define fi a (5m a ) as the probability distribution of the deviations of the local order 
parameter from its spatial mean 5m a (r) = m a (f) — (1/V) Jm a (f)df at p c a (L), for each 
a. A comparison of the averaged distributions \p, a ] for the two different lattice sizes (Fig. 
E]) reveals that the tail of the distribution, related to the significant deviations from the 
average (5m a ~ 1), becomes more pronounced with increasing system size, although the 
order parameter vanishes in the limit L — > oo. As shown in Fig. HJ the width of the 
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Figure 3: (color online) Distributions of spatial variations of the order parameter for two lattice sizes 
averaged over several disorder configurations (in parentheses) for c = 0.25. 
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Figure 4: (Color online) Size dependence of: a) order parameter \rn a ], b) width of the spatial variations of 
the local order parameter [|<5m a (r)|], and c) width of the thermal fluctuations of the order parameter [<5m a ], 
for c = 0.25. The dashed lines denote the best fits which include scaling corrections. 



spatial inhomogeneities of order parameter [|5m a (f)|] is larger than the width of thermal 
fluctuations of the order parameter [8m a ] = y[m 2 a — m^J. We conclude that spins which 
contribute to ordering are increasingly frozen as L increases. 

The magnetic exponent (-) calculated for c = 0.125 and 0.25 from the first moment of 
the order parameter (0.129(3)) which includes the contribution from inhomogeneities does 
not differ significantly from the value obtained from the second moment (0.132(3)), which 



includes only thermal fluctuations 



25j . Both of these values are within the range of previous 



results for the exponent f$/v 28rr30j]. We conclude that the difference in m a between different 



disorder configurations en, when determined at p c a , is of the same order as the width of the 



thermal fluctuations of the order parameter for each a. 

In conclusion, we have used the recently proposed EIC algorithm to study the critical 
behavior in the 2D q = 3 Potts model with quenched disorder by averaging at quasicritical 
temperatures of individual disorder configurations. This procedure has allowed us to sepa- 
rate the thermal from disorder-sample fluctuations and obtain the critical exponents y T and 
Dh = d — fi/v characterizing the transition in a single disorder configuration. The difference 
between the critical and disorder fluctuation exponents is interpreted by the emergence of 
spatial inhomogeneities on all scales at the transition point. By examining the local order 
parameter, we explicitly show that the inhomogeneities are present and that they become 
increasingly frozen (localized) as the size of the system increases. The above studies are 
currently being extended to higher dimensions, by using the same EIC algorithm. The ap- 
proach may be easily applied to more specific problems, such as the extensively discussed 2D 
Ising model [32| as a marginal case of quenched disorder according to the Harris criterion, 
or to the change from a first- to second-order phase transition induced by disorder [171 ]. 

Note added. Another interesting extension of the application of the EIC algorithm would 
include models with random field, for which an alternative approach to eliminate the lack 



of self-averaging has been recently proposed 33 ] . 
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